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Abstract 


This report presents the main results of the modeling task of the PPA project. The 
objective of this task is to make major progress towards developing a new computa- 
tional tool with new capabilities for simulating cylindrically symmetric 2.5 dimensional 
(2.5 D) PPA’s. This tool may be used for designing, optimizing, and understand- 
ing the operation of PPA’s and other pulsed power devices. The foundation for this 
task is the 2 D, cylindrically symmetric, magnetohydrodynamic (MHD) code PCAPPS 
(Princeton Code for Advanced Plasma Propulsion Simulation). PCAPPS was origi- 
nally developed by Sankaran (2001, 2005) to model Lithium Lorentz Force Accelerators 
(LLFA’s), which are electrode based devices, and are typically operated in continuous 
mode. The magnetic field in this code only has an azimuthal component. For the 
PPA project, PCAPPS was modified to apply to PPA’s, which are drive coil based 
devices operated in pulsed mode. The modification consists of changing the geometry 
of the computational domain, adding the radial component of the magnetic field to 
the model, and implementing a first principles, self-consistent algorithm to couple the 
plasma and power circuit that drives the plasma dynamics. This coupling algorithm 
is the most important result of the modeling effort so far. As in reality, the coupling 
occurs through the current and electric field in the drive coil. At a given time t the 
current in the drive coil provides a boundary condition for PCAPPS, allowing the code 
to be advanced to the next time step t + At. Then the voltage drop V p i a sma(t + At) 
across the drive coil due to the plasma electric field is computed and represented in the 
power circuit equations as an EMF in series with the drive coil’s resistance and induc- 
tance. The plasma electric field is computed from the Ohm’s law used in PCAPPS, 
and includes all effects of the plasma dynamics on the power circuit. The voltage drop 
across the drive coil is V p i asrna plus the voltage drops across the drive coil’s resistance 
and inductance. The power circuit equations are then solved to determine the drive coil 
current at time t + At, which provides an updated boundary condition for PCAPPS. 
This iterative algorithm is applied at each time step. The coupling algorithm appears 
to be new. It is general in that it allows any Ohm’s law, and any power circuit that 
couples to the plasma through a drive coil to be included in any PPA simulation. The 
basic structure of PCAPPS, including these modifications, is described. Simulations 
of pulsed inductive thrusters (PIT’s) having drive coil radii of 5.5 cm and 13.5 cm, and 
for Lithium and Argon propellants are presented. For these simulations, the thruster 
is powered by a capacitor in an LRC circuit using values of L,R, and C characteristic 
of TRW PIT thruster experiments. The objective of the modeling task was achieved 
in that: (1) An apparently new, self consistent, general, first principles based power 
circuit-plasma coupling algorithm was conceived and implemented; (2) Major modi- 
fications were successfully made to PCAPPS as part of the longer range objective to 
make the code fully 2.5 D. 



1.0 The Model 


PCAPPS is a two dimensional, cylindrically symmetric MHD code. Cylindrical coordinates 
(r, 9, z) are used, z is the coordinate along the axis of the cylinder, r is the distance from the 
cylinder axis, and all quantities are independent of 6. The drive coil is a disk in the plane 
z = 0. Its axis is the 2 axis, and the direction of increasing z is away from the drive coil. 
PCAPPS contains a description of the transport processes currently believed to be most 
important for the description of PPA’s. It contains the parallel (Spitzer) and Hall electrical 
conductivities, an expression for anomalous resistivity, scalar expressions for the electron 
and ion thermal conductivities, and a semi-empirical expression for radiative loss. 

The conservation equations for the total mass, momentum, and energy, and Faraday’s 
law are: 
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Here u and w are the radial and axial components of the velocity; Bg and B r are the 
azimuthal and radial components of the magnetic field; j r ,jz, and je are the components of 
the current density; p and p are the total mass density and total pressure; T e and T % are the 
electron and ion temperatures; n e and n % are the electron and total ion number densities; // 
is the sum of the Spitzer and anomalous resistivities; and e is the magnitude of the electron 
charge. Equations (4) and (5) are the nonzero components of Faraday’s law. The terms in 
these equations that are oc (en,) _1 are the Hall terms arising from the Hall conductivity 1 . 
K t h, e and Kth,i are the electron and ion thermal conductivities, //q is the magnetic permeability 

1 The exact Hall terms have n e in place of n t . The exact relation between the two is the quasi- neutrality 
condition n e = j n i> where the sum is over all ionization stages, and n{ is the number density of j 

times ionized ions. Then the approximation n e ~ n* = XlyLo n '\ is justified if the ionization occurs on time 
scales much shorter than those of interest, and if j is not much larger than unity, which appears to be the 
case. The reason for making the approximation is to avoid the very small computational time steps that 
result from n e < tij ~ nj at the beginning of the simulation before significant ionization occurs. One of the 



of vacuum. 7 is the ratio of specific heats for the gas of ions and electrons, and is given 
empirically as a function of T % in PCAPPS. The term n e riiL (T e ) describes radiative energy 

loss, where L(T e ) is an empirically determined function, e = \p{u 2 + w 2 ) + B$ ^ r + is 
the total energy density per unit volume. 7 is obtained from a model developed by Choueiri 
[3] that is valid for Lithium and Argon. The anomalous component of the total resistivity 
77, due to microinstabilities, is derived by Choueiri [4], The empirical function L(T e ) that 
occurs in the radiative loss term in the energy equation (6) is given by Post et al. [5] for a 
Lithium plasma. This function is also used in PCAPPS to describe radiative loss in Argon. 
It is found by Sankaran (2005) that radiative loss is not a significant cooling mechanism 
for LLFA plasmas. The same is true for the PPA plasmas considered here since they are 
characterized by similar temperature and density ranges. 

The terms involving B r and jo in equations (1) - (6) are not in the original version of 
PCAPPS. They are implemented here to extend the capability of PCAPPS to modeling 
PPA’s. As in the original version of PCAPPS, the 0 component of the momentum equation, 
and all terms involving B z are not included in the model. The extension of PCAPPS pre- 
sented here is the simplest one that allows for modeling PPA’s. An extension of PCAPPS 
to allow for modeling general 2.5 D, cylindrically symmetric plasma dynamics must include 
the 6 component of the momentum equation, and all terms involving B z . PCAPPS will be 
modified to include this equation and these terms in the second phase of the project during 
2006 - 2007. 

A simple analysis was done to estimate the error involved in omitting the 6 component 
of the momentum equation, and the terms involving B z when modeling PIT type devices. 
This analysis and its results are described as follows. The azimuthal component of the 
momentum equation may be written as dv/dt = ( j x B)e = Fg, where dv/dt is the convective 
(Lagrangian) time derivative of the azimuthal component v of the velocity. Similarly, the 
other components of the momentum equation may be written as du/dt = ( j x B ) r — dp/dr = 
F r ,dw/dt = ( j x B) 2 — dp/dz = F z . For a typical simulation of a PIT, described later in 
this report, it is found that \Fg\ ~ |F r | ~ 10 _3 |F 2 | at time t = 1.41 /i sec when the current in 
the power circuit reaches its first peak, and \F$\ ~ (10 -2 — 10 _1 )|F r | ~ \F Z \ at t = 7.25 p sec 
when the current in the power circuit reaches its second peak. This inequality is maintained 
for the remainder of the simulation. These results suggest that if the azimuthal component 
of the momentum equation is included in the model, only a relatively small azimuthal flow 
speed is generated, in which case the error in omitting this equation is small. The error in 
omitting terms containing B z is estimated by using the axial component of Faraday’s law, 
given by dB z /dt = — (V x E) 2 , with E given by the Ohm’s law of the model. For the 
same PIT simulation, V x E is computed, and its r, 6, and z components are compared. 
The result is that |(V x E) z | ~ 10 -1 |(V x E) o\ r '- / | ( V x E) r | at t — 1.41 p sec, and 
|(V x E) 2 | ~ |(V x E)o| ~ (lO -2 - 10 -1 )|(V x E) r | at t — 7.25 /i sec. These results suggest 
that for a PIT type device the effects of B z are small compared with those of B r , but may 
be comparable with those of Bq. This is consistent with experimental results for the TRW 
PIT showing that almost no axial magnetic field is generated during device operation [6, 7]. 


tasks of the follow on phase of the PPA project is to explore the error introduced by this approximation, and 
compare PCAPPS and MACH2 with respect to how they treat the Hall terms, which can have important 
effects in 2 and 3 dimensions. 
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Equations (l)-(6) must be supplemented by additional equations to complete the model. 
These equations are the electron energy equation, the ideal gas equations of state for the 
electrons and ions, the caloric equation of state for the electrons, and Saha’s equation relating 
the number densities of electrons and the several ionization stages of the ions. 

The electron energy equation is given by 
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Here £ e ,p e ,m e , Mi,k B , and v ei are the electron energy density per unit volume, electron 
pressure, electron mass, ion mass, Boltzmann’s constant, and the electron-ion collision fre- 
quency. The first, second, and third terms on the right hand side of equation (7) correspond 
to Ohmic heating, electron-ion energy exchange via collisions, and electron heat flux. 

The ideal gas and caloric equations of state for the electrons are 
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Here y e = 5/3. 

Saha’s equation is 


wjw e _ 2 (27rm e ft; B r e ) 3/2 Eiff/e 
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Here j and j — 1 label ionization stages, l labels discrete atomic energy levels, and gi is the 
statistical weight of the corresponding energy level for a given ionization stage. The use of 
this Saha equation is based on the assumptions that the distribution of the ions over their 
internal energy levels (bound and continuum (ionized)) is determined entirely by electron- 
ion collisions, that this distribution is an equilibrium (Boltzmann) distribution, that the 
electrons are in equilibrium at temperature T e , and that the distribution of the ions over 
their internal energy levels is in equilibrium with the electrons. The electrons and ions are 
assumed to respond instantaneously on MHD and radiative time scales to maintain these 
equilibrium conditions. A discussion of the conditions under which these assumptions are 
valid is given in Sankaran (2001,2005) and references therein. Roughly, the assumptions are 
valid if r ee <C T m -C MHD and radiative transition time scales. Here r ee and r e , are the 
electron-electron and electron-ion collision times. 

The ion pressure, and equation of state are 

Pi=P~Pe (10) 

Ti = ( 11 ) 

riiKs 

Complete details of the original PCAPPS model and code, including a discussion of the 
various approximations, are given by Sankaran (2001, 2005). 


3 



1.1 Power Circuit - Plasma Coupling Algorithm 

One of the most challenging problems in modeling any inductive electric thruster is to model 
the coupling between the power circuit and plasma with sufficient accuracy. This problem 
demands a self consistent treatment for an accurate solution since the power circuit-plasma 
coupling is strongly bi-directional. The inductive electric field generated by the power circuit 
ionizes and accelerates the propellant. The resulting plasma dynamics generates its own 
electric field that reacts back on the power circuit, effectively having inductive and resistive 
effects on the circuit. 

Thus far, there have been two basic approaches to modeling the power circuit - plasma 
coupling. One approach is to treat the plasma as an LRC circuit inductively coupled to the 
power circuit [6], [7]. This approach includes bi-directional coupling, but the effect of the 
complex plasma dynamics on the power circuit is treated approximately in terms of lumped 
circuit elements. The second approach is to drive the plasma with an experimentally deter- 
mined drive coil current waveform [8]. Mikellides and Neilly perform MACH2 simulations 
of a PIT driven by such a waveform. This waveform is used to determine L, R, and C in 
the equivalent LRC circuit that reproduces the current waveform. This approach, though 
based on measurements of the current in an operating PIT, is uni-directional, and includes 
all effects of the power circuit and the plasma dynamics together through the values of L, R, 
and C . These restrictions of uni-directional coupling and lumped element circuit modeling 
of the effects of plasma dynamics are significant, indicating the need for a bi-directional 
coupling algorithm that includes the exact plasma dynamics. 

Here a simple, self consistent, first principles based algorithm for including the bi- 
directional power circuit-plasma coupling in PPA simulations is presented and implemented. 
The algorithm is as follows. The effect of plasma dynamics on the power circuit must occur 
through the plasma electric field E acting on the charges in the drive coil. The effect of E is 
to alter the voltage Vd between the terminals of the drive coil from what it is in the absence 
of the plasma. The contribution V p i asma of E to Vi is exactly computed as the scalar line 
integral of E along the wires between the terminals of the drive coil. V p i asma is then included 
as an EMF in series with the drive coil resistance and inductance in the circuit equations 
governing the time evolution of the current through the drive coil in the power circuit. Given 
Vpiasma at time t, the current I(t + At) in the drive coil is computed by solving these equa- 
tions. I{t + At) is then used to determine an updated current density at the surface of 
the drive coil, which provides an updated boundary condition for PCAPPS. This iterative 
procedure is performed at each time step. E contains all effects of the plasma dynamics 
on the power circuit, including time dependent resistive and inductive effects, regardless 
of the complexity of the plasma dynamics. Vd equals the sum of V p i asma and the voltage 
drops across the R and L in the power circuit that represent the drive coil’s resistance and 
inductance, which may be measured, and are independent of time. 

For an arbitrary power circuit the algorithm is easily implemented using a circuit mod- 
eling software package such as PSPICE. The drive coil is represented in any power circuit 
as a constant inductance and resistance, and a time dependent V p i asma , all in series. The 
algorithm iterates between PSPICE and PCAPPS as described above, with PSPICE provid- 
ing an updated drive coil current boundary condition for PCAPPS, and PCAPPS providing 
Vpiasma at each time for PSPICE to compute a new drive coil current. 
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Additional insight into the coupling algorithm is gained by using the electromagnetic 
boundary condition which states that the tangential component of the total electric field is 
continuous across any material interface in its rest frame. Let E unre be the electric field that 
drives and is parallel to the current in the drive coil. Let E p be the component of the plasma 
electric field at the surface of, and parallel to the drive coil wires. The boundary condition 
requires that E vnre = E p + E' , where E' is the electric field at the boundary between the wires 
and the plasma volume that contains the effects of the drive coil’s constant resistance and 
inductance on Vg. The scalar line integral of E' along the wires between the terminals of the 
drive coil equals the voltage drop across the drive coil due to this resistance and inductance. 

For the simulations presented in this report the power circuit is an LRC circuit. The 
values of L, R, and 6', and other simulation parameters are presented in §3.2. Fig. 1 is a 
schematic diagram of the power circuit-drive coil configuration. The drive coil is shown as 
a set of concentric wire rings. The axis of the rings is the z axis of the simulation, the 
drive coil occupies the area defined by (z = 0, 0 < r < R co u), and the wires run along the 
0 direction. This approximation derives from the design of the original TRW PIT drive 
coil, which is comprised of 18 spiral strands of wire that can be accurately approximated 
as 9 concentric wire rings, each carrying an azimuthal current /. For the simulations, the 
drive coil is assumed to be a continuous disk carrying a uniform (r independent), azimuthal 
current density jg given by 

jo = (12) 

where A is the part of the cross sectional area of each wire in the drive coil through which 
current flows. It is given by 

A = tiR 2 w - it (R w - S) 2 , (13) 

where R w is a characteristic radius of a wire strand (R w ~ 1 cm), and S is the skin depth 
of copper at the natural frequency of the LRC circuit uncoupled from the plasma. Then 
the current density at each of the guard cells in the simulation domain along the drive coil 
radius 0 < r < R cm i at z = 0 is assigned the value jg. jg is updated at each time step, and 
provides the MHD boundary condition at the drive coil that allows PCAPPS to compute 
the state of the plasma at the next time step. 

The circuit equations for the power circuit - plasma system are: 


dl 

(14) 

V = RI + L— + V p i asma 

Q = CV 

(15) 

dQ _ j 
dt 

(16) 


where V is the voltage across the capacitor, and V p i asma is obtained from: 


,0 _ _ N 

Vplasmaij) = E • dl — 27 T Y,RiE e (R u t). 

J Router A — 1 


Here Eg is the component of the plasma electric field at, and parallel to the drive coil wire 
rings, Ri is the radius of guard cell i from the z axis, and N is the number of guard cells 
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Conceptual Circuit Plasma Coupled System 


-Wv 


MJiJ 



Figure 1: Schematic of the Power Circuit - Plasma System 


along the radial direction. Eg is computed using the PCAPPS Ohm’s law as 

Eg = ( -v x B + 7 ] j + - - B ) . (18) 

V ne Jo 

The terms on the right hand side of this equation are the convective, resistive, and Hall 
contributions to Eg. They are updated at each time step. 

The capacitor voltage V, and circuit current / are updated by solving equations (14) - 
(16). The solution to these equations is 


V(t) 


nt) 


+ 


+ 

+ 


r x % 

„—\t 


A 


c o 


sin (cut) + cos (cut) 


sin (cut) 


e Xt ' ( Vpiasma (t*)) 


ljLC 


e xt cos (c ot) ft xt 


uLC 


e M ' ( Vpiasma (f)) 


cos (uit 1 ) dt' 
sin (c ot') dt' 


CV t ,e 


-xt 


A 2 + io 2 . 


LU 


sin (cot) 


e xt (— A sin (ut) + co cos (cot)) 

loL J 

e~ xt (A cos (ut) + lo sin (utt)) r* 
ojL Jo 


f e M ' (Vpiasma (t'))cOs(u>t')dt' 
e x, ‘ (Vpiasma ((')) sin (wt') dt ' , 


(19) 


( 20 ) 


where 

= a “ dA = ^' (21) 
I determines the current density at the drive coil that is used as a boundary condition for 
PCAPPS. The initial conditions 1(0 ) = 0, P(0) = Vq are assumed. PCAPPS is initialized 
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for its first time step by computing dl /dt at t = 0, and multiplying the result by a seed time 
step A t s to obtain the required drive coil current boundary condition. For the simulations 
presented in §3 the choice A t s = 5 x 10 -9 sec is made. The Courant condition then modifies 
this time step, if necessary, and determines all subsequent time steps. A flow chart of the 
power circuit - plasma coupling algorithm is shown in Fig. 2. 


2.0 Numerical Algorithm 

The model uses a characteristic splitting finite volume method. In this method, in each cell 
j, each fluid equation reduces to 


iu--wJ F ' dl+<s>i - (22) 

Here fl, is the cell volume (area in 2D), the line integral is around the boundary of cell j, 
Uj is a typical conserved variable (e.g. p), F is the corresponding flux, and < S >j is the 
average of the source term, when present, over fl,. 

There exists a large literature on finite volume methods. For a comprehensive, thorough 
treatment see Godlewski & Raviart [9]. These methods have emerged as powerful methods 
in hydrodynamics and MHD since they employ semi-analytic techniques that allow for the 
resolution of strong discontinuities without the use of filtering schemes such as flux-corrected 
transport (FCT). Equation (22) is obtained by multiplying each fluid equation by a trial 
function, and integrating over all space. Here, and in most cases, the trial function is chosen 
to have unit value in a given cell j, and be zero otherwise. Then the integration over all space 
reduces the calculations to those involving a single cell j and its boundaries, as indicated in 
equation (22). Analytic calculations are then employed to determine the values of Uj at cell 
boundaries in terms of its values at cell centers using a method developed by Roe [10] . This 
method is based on the analytic Rankine-Hugoniot jump conditions. 

The following example shows the effectiveness of the Roe method in resolving strong 
discontinuities, and serves as a test case for validating the modification of PCAPPS. This 
example consists of triggering a fast MHD shock wave by using a step function in B r as 
an initial condition, and then checking the degree to which the Rankine-Hugoniot jump 
conditions across the shock front are satisfied. This step function is a change in B r from 500 
to 100 G at the drive coil. It models pulsing a PIT type device with a very short current rise 
time. A shock wave develops as expected, with all profiles steepening into discontinuities 
that travel together at the same speed. The Rankine Hugoniot jump condition for each 
profile is checked. For B r this condition is obtained from the Faraday law equation (5), and 
is given by 

sAB r = A F z . (23) 

Here s is the shock speed, determined from the speed of the discontinuity in the simulation, 
and A B r and A F z are the change of B r and its flux across the discontinuity. The Rankine- 
Hugoniot conditions are found to be satisfied to within an error ~ 15%. This error is 
considered small since the simulation is 2 D while the Rankine-Hugoniot jump conditions 
are derived for a 1 D system. 
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Flow Chart of Circuit Coupling Process 
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Figure 2: Flow Chart of the Power Circuit - Plasma Coupling Algorithm. 
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An additional test of the code modification of adding B r to the model is to check the 
coupling of B r and Bq by the Hall electrical conductivity. The test consisted of initializing 
B r as a sine wave, and then checking that a corresponding Bq is generated as a sine wave. 
This was found to be the case. The two sine waves have approximately the same amplitude 
in the absence of the plasma resistivity 77, as expected. 

3.0 Examples of Simulations of Pulsed Plasma Acceler- 
ators 

3.1 Geometry 

The original version of PCAPPS uses a mesh of quadrilaterals to fit an LLFA geometry based 
on the device design by Polk [11]. Here the geometry is modified to that of a PIT. The radius 
Rcoii °f the drive coil is chosen to be 5.5 cm or 13.75 cm. The first value corresponds to a 
drive coil area about 1.7 times smaller than that of a FARAD thruster [12]. The second value 
corresponds to a drive coil area about 11 times smaller than that of the TRW PIT drive coil, 
which is an annulus with an inner radius of ~ 20 cm, and an outer radius of ~ 50 cm [6, 7] . 
Therefore, the simulations presented here are for compact devices. It is important to note 
that Rani is not the radius of the inlet, over which the propellant is uniformly distributed. 
The inlet radii corresponding to the smaller and larger values of R coi i are R, n i e t = 11 cm 
and 27.5 cm, respectively. Simulations using TRW PIT size drive coils that cover the entire 
inlet area are currently being performed 2 . As expected, they show much stronger coupling 
between the power circuit and plasma, leading to higher I sp . Fig. 3 shows the geometry and 
its mesh. The figure shows an r — z plane cut through the cylindrically symmetric device. 
The 2 axis is the axis of the device. The disk shaped drive coil is at the left end of the 
device, normal to the page. The 2 axis passes through the center of the drive coil at z = 0. 
The current density at the drive coil is used as a boundary condition for the Faraday law 
equations for B r and Bq, allowing PCAPPS to advance the magnetic field in time. All other 
boundaries are treated as free stream boundaries. 

3.2 Simulation Parameters 

The values L = 100 nH, C = 9 77/, R = 6.7 x 10 -3 O, and bo = 30 kV are chosen for all 
simulations. These values are characteristic of TRW PIT experiments [7]. It is assumed the 
propellant fills the entire computational volume with an initial density n 0 = 10 21 m -3 . This 
value may be compared with the value Uq ~ 8.7 x 10 21 m -3 corresponding to actual TRW 
PIT experiments in which a 3 cm thick layer of 2 mg of Li propellant is distributed over 
a drive coil area ~ 0.66 m 2 . The initial ionization level is assumed to be 5%. Simulations 
are performed using Li and Ar propellants, the difference between the two cases being the 
corresponding equations of state and mass densities. The assumption of propellant filling 
the entire computational volume is not realistic. Simulations using a more realistic initial 

2 The inlet is a remnant of the original LLFA geometry in PCAPPS. Current simulations, not described 
in this report, effectively replace the inlet by the drive coil, and extend the radius of the drive coil to the 
boundary of the computational domain. 
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Design to emulate TRW/PIT &/or FARAD 

— interface of the plasma and solid boundaries 
Centerline (cylinder axis) 

— No radial convective fluxes, gradients, or thermal conduction 
Inlet: Drive coil is tangent 

— Where the mass flow of the propellant takes place 
Free Stream 

— No normal gradients 

— Material flows in and out freely 

— All other boundaries are free stream boundaries 
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Figure 3: Device geometry, r and z are in meters. 


density profile that decreases exponentially fast with increasing distance from the drive coil 
are currently being performed. 

3.3 Simulation 1: R co u = 5.5 cm, Uni-Directional Power Circuit- 
Plasma Coupling, Li Propellant 

Uni-directional coupling means the power circuit acts on the plasma, but does not receive 
any feedback from the plasma. This coupling is defined by setting V p i asma = 0 in the 
power circuit equations. This simulation is performed to allow comparison of results for uni- 
directional and bi-directional coupling. As expected, there are major differences, emphasizing 
the importance of using a bi-directional power circuit-plasma coupling algorithm. 

Figures 4 and 5 show the circuit current I dr it) and plasma voltage V p i asma (t) for a run 
lasting 50 /i sec. I cir shows oscillations with a period ~ 5.9 /j sec, which is the natural 
frequency u of the LRC circuit. For the chosen values of L, R, and C, oj 

The magnitude of the current rises from zero to 2.65 x 10 5 A in the first quarter cycle 
as expected. It decreases by ~ 18% between successive peaks due to the circuit’s resistance. 
B r reaches a maximum value of 250 G at the drive coil. As shown by the next example in 
§3.4, this dissipation rate substantially increases when the power circuit-plasma coupling is 
bi-directional. V p i asma is periodic at approximately the same frequency as the power circuit. 
It rises to ~ 2200 V which is more than one order of magnitude below the initial capacitor 
voltage of 30 kV. This relatively small value of V p i asma is attributed to the relatively small 
power circuit-plasma coupling area of 9.5 x 10 -3 m 2 , which is the area of the drive coil, 
compared to a coupling area 0.66 m 2 for the TRW PIT drive coil. 

Figures 6, 7, and 8 show the electron density at three successive times. The initial blue 
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Figure 4: Circuit current (A) vs. t(sec); R co u = 5.5 cm, Uni-directional Coupling, Li 
Propellant. 
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Figure 5: V p i asma (V) vs. t (sec); R co u = 5.5 cm, Uni-directional Coupling, Li Propellant. 
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background corresponds to the uniform initial background electron density of 5 x 10 19 m 3 , 
which is the initial ionization fraction of 5 % times n 0 • Figures 9, 10 and 11 show the 
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Figure 6: Electron number density (m 3 ) at t — 6.9// seconds; R co u = 5.5 cm, Uni-directional 
Coupling, Li Propellant; r-z in meters. 



Figure 7: Electron number density (rri 3 ) at t = 25 // sec; R co n = 5.5 cm, Uni-directional 
Coupling, Li Propellant; r-z in meters. 

total number density at the same times. Figures 12, 13, and 14 show the corresponding 

axial velocities (here labeled V z ), and figures 15, 16, and 17 show the corresponding radial 
velocities (here labeled V r ). The large value of V z near the inlet at later times is a result 
of the power circuit continuing to drive the plasma 3 . The profiles show plasma compression 
and acceleration from the drive coil to the exit on the far right. The maximum I sp for this 
simulation is 3400 sec. 

3 Subsequent simulations show that using a sufficiently large drive coil radius causes almost all of the 
energy of the power circuit to be transferred to the plasma in approximately one oscillation of the circuit 
current. 
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Figure 8: Electron number density (rri 3 ) at t = 45 // sec; R co ,i = 5.5 cm, Uni-directional 
Coupling, Li Propellant; r-z in meters. 



Figure 9: Total number density ( m 3 ) at t = 6.9 // seconds ]R CO ii = 5.5 cm, Uni-directional 
Coupling, Li Propellant; r-z in meters. 


13 



Figure 10: Total number density (rri 3 ) at t = 25 fi seconds; R co u = 5.5 cm, Uni-directional 
Coupling, Li Propellant; r-z in meters. 
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Figure 11: Total number density ( m 3 ) at t = 45 /.( seconds; B co h = 5.5 cm, Uni-directional 
Coupling, Li Propellant; r-z in meters. 
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Figure 12: Axial velocity (m — sec x ) at t = 6.9 fi seconds ]R CO ii = 5.5 cm, Uni-directional 
Coupling, Li Propellant; r-z in meters. 
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Figure 13: Axial velocity ( m — sec x ) at t = 25 fi seconds; R co u = 5.5 cm, Uni-directional 
Coupling, Li Propellant; r-z in meters. 
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Figure 14: Axial velocity (m — sec x ) at t = 45 /r seconds ;Rcoii = 5.5 cm, Uni-directional 
Coupling, Li Propellant; r-z in meters. 



Figure 15: Radial velocity (m — sec x ) at 6.9 /i seconds ]R co u = 5.5 cm, Uni-directional 
Coupling, Li Propellant; r-z in meters. 
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Figure 16: Radial velocity (m — sec x ) at t = 25 /r seconds ]R co u = 5.5 cm, Uni-directional 
Coupling, Li Propellant; r-z in meters. 
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Figure 17: Radial velocity (m — sec x ) at t = 45 fi seconds; R CO j/ = 5.5 cm, Uni-directional 
Coupling, Li Propellant; r-z in meters. 
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3.4 Simulation 2: R co u = 5.5 cm, Bi-Directional Power Circuit- 
Plasma Coupling, Li Propellant 

This simulation differs from Simulation 1 only in that there is now bi-directional coupling 
between the power circuit and the plasma. Figures 18 and 19 show I C i r (t) and V p i asma (t) for 
a period of 50 /j, seconds. A comparison of Fig. 4 and Fig. 18 shows a significantly larger 
damping rate for the case of bi-directional coupling. The current in Fig. 4 decreases by 
~ 18% per period. The current in Fig. 18 decreases by 30 — 40% during each successive 
oscillation. The oscillation frequency for I cir is essentially equal to the natural frequency u 
of the power circuit. Since ui ~ (L(7) _1//2 for the chosen values of L, R, and (7, this implies 
that the inductive coupling between the power circuit and the plasma is not significant. 
In reality, this coupling can be significant. It is conjectured that the lack of inductive 
coupling in this simulation, and in the other simulations discussed below, is mainly due to 
two factors. The first is the unrealistic feature of the simulation that the propellant fills 
the entire computational volume, instead of the initial propellant distribution being a thin 
layer of gas on the drive coil. The second factor is that the mutual inductance between the 
plasma and the power circuit is expected to be proportional to the coupling area, which is 
the area of the drive coil. Then the coupling is oc R 2 d . The maximum value of R,i used in the 
simulations is 13.75 cm. It may be that this value is not large enough to cause significant 
inductive coupling. Simulations are currently being performed to resolve this issue. 



Figure 18: Power circuit current (A) vs. t (sec). R co u = 5.5 cm, Bi-directional Coupling, Li 
Propellant. 

Figures 20, 21, and 22 show the electron density at three successive times. 

Figures 23, 24, and 25 show the total number density at the same times. Figures 26, 27, 
and 28 show the corresponding axial velocities, and figures 29, 30, and 31 show the corre- 
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Figure 19: V p i asma (V) vs. t (sec). R co u = 5.5 cm, Bi-directional Coupling, Li Propellant. 



Figure 20: Electron number density (m 3 ) at t = 6.9/U seconds. R co u = 5.5 cm, Bi-directional 
Coupling, Li Propellant; r-z in meters. 
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Figure 21: Electron number density (m 3 ) at t = 25 fi seconds. R co u = 5.5 cm, Bi-directional 
Coupling, Li Propellant; r-z in meters. 



Figure 22: Electron number density (m 3 ) at t = 45 fi seconds. R co u = 5.5 cm, Bi-directional 
Coupling, Li Propellant; r-z in meters. 
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sponding radial velocities. The maximum I sp achieved for this simulation is 2500 seconds. 
This value is significantly smaller than the corresponding value of 3400 seconds for Simula- 
tion 1. This reduction in the maximum value of I sp , and the increase in the dissipation rate 
of the current in the power circuit over its value for Simulation 1 is believed to be due to the 
fact that, for bi-directional coupling, electrical energy in the plasma can be dissipated by the 
resistance of the power circuit, as well as by the plasma’s own resistance. When the coupling 
is uni-directional, plasma electrical energy can only be dissipated by the plasma’s resistance. 
Bi-directional coupling opens up a new channel for the resistive dissipation of plasma elec- 
trical energy, resulting in more energy being thermalized, and less being converted into bulk 
flow kinetic energy, leading to a lower I sp . 
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Figure 23: Total number density (rn 3 ) at t = 6.9 /.( seconds; R co n = 5.5 cm, Bi-directional 
Coupling; r-z in meters. 



Figure 24: Total number density (m 3 ) at t = 25 /j seconds; R co u = 5.5 cm, Bi-directional 
Coupling, Li Propellant; r-z in meters. 
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Figure 25: Total number density (m 3 ) at t = 45 /j seconds; R co n = 5.5 cm, Bi-directional 
Coupling, Li Propellant; r-z in meters. 



Figure 26: Axial velocity ( m — sec 3 ) at t = 6.9 seconds; R co u = 5.5 cm, Bi-directional 
Coupling, Li Propellant; r-z in meters. 
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Figure 27: Axial velocity (m — sec x ) at t = 25 n seconds; R co u = 5.5 cm, Bi-directional 
Coupling, Li Propellant; r-z in meters. 


V z: -lOOOO -1200 7600 16400 25200 


0.1 

O 



Figure 28: Axial velocity (m — sec x ) at t = 45 /r seconds; R co n = 5.5 cm, Bi-directional 
Coupling, Li Propellant; r-z in meters. 
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Figure 29: Radial velocity (m — sec x ) at t = 6.9 /i seconds; R co u = 5.5 cm, Bi-directional 
Coupling, Li Propellant; r-z in meters. 



Figure 30: Radial velocity ( m — sec x ) at t = 25 fi seconds; R co u = 5.5 cm, Bi-directional 
Coupling, Li Propellant; r-z in meters. 
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Figure 31: Radial velocity (rri — sec x ) at t = 45 fi seconds; R co u = 5.5 cm, Bi-directional 
Coupling, Li Propellant; r-z in meters. 

3.5 Simulation 3: R co u = 13.75 cm, Bi-Directional Power Circuit- 
Plasma Coupling, Li Propellant 

Here the drive coil-plasma coupling area is increased by a factor of 6.25 over that for Sim- 
ulation 2. Idr(t ) and V p i osma (t) are shown in figures 32 and 33. A comparison of figure 
32 with figure 18 shows that the current damps at about the same rate as for Simulation 2. 
This is believed to be mainly due to the result of the unrealistic assumption of a constant 
drive coil resistance R. Since the area of the drive coil is increased by the factor 6.25 over 
that in Simulation 2, the resistance should increase by approximately the same factor. This 
is expected to increase the current damping rate by about the same factor. A comparison 
of figure 33 with figure 19 shows that V p i asma decreases more rapidly than for Simulation 2. 
This is believed to be due to the larger coupling area allowing for more rapid dissipation of 
plasma electrical energy by the power circuit resistance. 

Figures 34, 35, and 36 show the total density at three successive times. 
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Figure 32: Circuit current (A) vs. t (sec). R co u = 13.75 cm, Bi-directional Coupling, Li 
Propellant. 
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Figure 33: Plasma voltage (V) vs. t (sec). R co u = 13.75 cm, Bi-directional Coupling, Li 
Propellant. 


26 





Figure 34: Total number density (m 3 ) at t = 10/U seconds. R co u = 13.75 cm, Bi-directional 
Coupling, Li Propellant; r-z in meters. 
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Figure 35: Total number density (m 3 ) at 7 = 25 // seconds. R co u = 13.75 cm, Bi-directional 
Coupling, Li Propellant; r-z in meters. 


Figures 37, 38, and 39 show the electron density at the same times. 

Figures 40, 41, and 42 show the corresponding axial velocity. 

Figures 43, 44, and 45 show the corresponding radial velocity. The duration of this 

simulation is 75 /i seconds. The increase in runtime over that of Simulation 2 is necessary to 
fully resolve the plasma dynamics since the axial velocity is much smaller than for Simulation 
2. This decrease in axial flow velocity is believed to be due to the greater dissipation rate 
of plasma electrical energy by the resistance R of the power circuit, thereby decreasing 
the energy available for propellant acceleration. This is a consequence of the bi-directional 
coupling of the power circuit and plasma allowing plasma electrical energy to be dissipated 
by the circuit resistance as well as by the plasma resistance. The maximum value of I sp for 
this simulation is 620 seconds. 


27 




Figure 36: Total number density (m 3 ) at t = 45 // seconds. R co ,i = 13.75 cm, Bi-directional 
Coupling, Li Propellant; r-z in meters. 



Figure 37: Electron number density ( m 3 ) at t = 10 /r seconds. R co u = 13.75 cm, Bi- 
directional Coupling, Li Propellant; r-z in meters. 
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Figure 38: Electron number density (m 3 ) at t = 25 p seconds. R co u = 13.75cm, Bi- 
directional Coupling, Li Propellant; r-z in meters. 



Figure 39: Electron number density (m 3 ) at t = 45 p seconds. R co u = 13.75cm, Bi- 
directional Coupling, Li Propellant; r-z in meters. 
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Figure 40: Axial velocity (m — sec x ) at t = 10 fi seconds. R co u = 13.75 cm, Bi-directional 
Coupling, Li Propellant; r-z in meters. 



Figure 41: Axial velocity (m — sec 1 ) at t = 25 /.i seconds. R co u = 13.75 cm, Bi-directional 
Coupling, Li Propellant; r-z in meters. 
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Figure 42: Axial velocity (m — sec x ) at t = 45 fi seconds. R co u = 13.75 cm, Bi-directional 
Coupling, Li Propellant; r-z in meters. 



Figure 43: Radial velocity (m — sec x ) at t = 10 fi seconds. R co u = 13.75 cm, Bi-directional 
Coupling, Li Propellant; r-z in meters. 
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Figure 44: Radial velocity (rri — sec x ) at t = 25 fi seconds. R co u = 13.75 cm, Bi-directional 
Coupling, Li Propellant; r-z in meters. 



Figure 45: Radial velocity ( m — sec x ) at t = 45 n seconds. R co u = 13.75 cm, Bi-directional 
Coupling, Li Propellant; r-z in meters. 
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3.6 Simulation 4: R co u = 13.75 cm, Bi-Directional Power Circuit- 
Plasma Coupling, Ar Propellant 

This simulation is identical to Simulation 3 except that Ar is used as a propellant instead of 
Li. The masses of Ar and Li are 39.948 amu and 6.94 amu. The difference between the Li 
and Ar cases enters the model through the mass density, and equation of state. 

Figures 46 and 47 show I dr it) an d V p i asma (f). Comparison of these figures with figures 32 
and 33 for Simulation 3 show that all profiles oscillate with a period essentially equal to 
the natural frequency u ~ 5.9// sec of the power circuit, but that the amplitude of I cir is 
significantly larger for the case of Ar. As for Simulation 2, the fact that the oscillation 
frequency is u, which is ~ (LC) -1 / 2 , indicates there is no significant inductive coupling 
between the plasma and the power circuit for these simulations. The result that the damping 
rate for I ar is smaller in the case of Ar than in the case of Li is due to the fact that 
relatively little work is done by the Lorentz force on the plasma during the simulation with 
Ar compared with that done during the simulation with Li. For a given axial Lorentz force 
F z , the acceleration of Ar is less than that of Li by roughly the mass ratio. The power 
applied by the power circuit to accelerate the plasma is ~ F Z V Z , which is smaller for Ar than 
for Li due to the smaller value of V z . 

For given F z , the ratio of the maximum values of I sp for the Li and Ar cases should be 
approximately equal to the inverse ratio of their masses. The maximum value of I sp for this 
simulation, obtained from figure 50, is 140 seconds. Then these ratios are 4.43 and 5.64. 

As a result of the increase in mass of the propellant, the plasma only moves a distance 
~ 0.2 m during the 75 /i sec duration of the simulation. 



Figure 46: Circuit current (A) vs. t (sec). R co u = 13.75 cm, Bi-directional Coupling, Ar 
Propellant 
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Figure 47: Plasma voltage (V) vs. t (sec). R co u = 13.75 cm, Bi-directional Coupling, Ar 
Propellant. 

Figures 48, 49, 50, and 51 show the total number density, electron number density, 
and axial and radial velocities at t = 75 /r sec. Increasing the initial capacitor voltage Vo 
increases the acceleration. However, ongoing work shows that using a realistic initial density 
profile, a significantly higher I ap is obtained for Ar without changing Vq. 
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Figure 48: Total number density (m 3 ) at t = 75 // seconds. R co n = 13.75 cm, Bi-directional 
Coupling, Ar Propellant; r-z in meters. 


1 .7581 4E+1 9 



O 0.2 0.4 0.6 


Figure 49: Electron number density ( m 3 ) at t = 75 fi seconds. R co u = 13.75 cm, Bi- 
directional Coupling, Ar Propellant; r-z in meters. 


4.0 Summary 

PCAPPS was modified to model pulsed inductive thrusters by: (1) Changing the geometry 
of the computational domain from that corresponding to an LLFA to that corresponding to a 
PPA; (2) Adding the radial component of the magnetic field to the model; (3) Implementing a 
general, first principles based power circuit-plasma coupling algorithm. These modifications 
are the first in a series of modifications to generalize PCAPPS to a complete 2.5 D simulation 
tool for cylindrically symmetric PPA’s, and similar pulsed power devices. 

An accurate algorithm for modeling power circuit-plasma coupling in PPA simulations 
is necessary for the simulation to be useful for designing, optimizing, and predicting the op- 
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Figure 50: Axial velocity (m — sec x ) at t = 75 fi seconds. R co u = 13.75 cm, Bi-directional 
Coupling, Ar Propellant; r-z in meters. 



Figure 51: Radial velocity (rri — sec x ) at t = 75 n seconds. R co u = 13.75 cm, Bi-directional 
Coupling, Ar Propellant; r-z in meters. 
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eration of PPA’s. Previous power circuit-plasma coupling algorithms are uni-directional or 
represent the plasma as a circuit coupled to the power circuit through a mutual inductance. 
For the PPA project, ISR conceived and implemented a self consistent, general, first prin- 
ciples based power circuit-plasma coupling algorithm. This algorithm appears to be new. 
The implementation was done using an L, R, C power circuit, and the Ohm’s law in the 
original version of PCAPPS. However, the algorithm is easily implemented for an arbitrary 
power circuit and Ohm’s law. For a power circuit that is not easily represented as a set 
of differential equations that may be solved analytically or numerically, a general circuit 
modeling tool such as PSPICE may be used to solve for the drive coil current that provides 
the required boundary condition for PCAPPS. PCAPPS then provides an updated voltage 
across the drive coil due to the plasma electric field. This voltage drop is represented as an 
EMF in series with the drive coil resistance and inductance in the power circuit, allowing 
PSPICE to compute an updated drive coil current. Changing the Ohm’s law in PCAPPS 
simply changes the algebraic equation for the plasma electric field. 

Since the end of the PPA project reported on here, PCAPPS was further modified in 
two important ways. The first was to expand the drive coil area to cover the entire area 
of the computational boundary at the back wall z = 0. This provides a more accurate 
representation of the geometry of real PPA’s The second modification was to implement 
a realistic initial propellant density profile, modeled as a density distribution that falls off 
exponentially with increasing distance from the drive coil. Simulation results incorporating 
these modifications will be reported during the current phase of the project in 2006 - 2007. 
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